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A numerical model to predict chill down in cryogenic 
transfer lines has been developed. Three chill down cases 
using hydrogen as the working fluid are solved: 1) a simplified 
model amenable to analytical solution, 2) a realistic model of 
superheated vapor flow, and 3) a realistic model of initially 
subcooled liquid flow. The first case compares a numerical 
model with an analytical solution with very good agreement 
between the two. Additionally, the analytical solution provides 
a convenient way to look at parametric effects on the chill 
down. The second and third cases are numerical models which 
provide temperature histories of the fluid and solid tube wall 
during chill down as well as several other quantities of interest 
such as pressure and mass flow rate. Of great interest is the 
ability to predict accurate values of chill down time (the time 
required to achieve steady-state cryogenic flow). The models 
predict that a 26 in. long, 3/16 in. ID aluminum tube has a 
shorter chill down time (-100 sec) and uses less hydrogen with 
superheated vapor flow than with initially subcooled liquid 
flow (>200 sec for chill down). 
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= heat transfer area, ft 
= tube cross-sectional area, ft 
= specific heat of the fluid, Btu/lb m -°R 
= specific heat of the tube wall, Btu/lb m -°R 
= specific heat at constant pressure, Btu/lb m -°R 
= tube inner diameter, ft 
= advected energy, Btu 
= thermally conducted energy, Btu 
= Darcy friction factor 
= mass flux, lb m /ft 2 -hr 
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= gravitational constant, 32. 17 ft-lbm/lbrs" 

= heat transfer coefficient, Btu/hr-fr-°R 
= modified Bessel function 
= specific enthalpy, Btu/lb m 

= conversion factor for mechanical work to heat, 778. 16 ft-lbj/Btu 
= flow resistance coefficient, lbps /lb„f -fr 
= vapor phase thermal conductivity, Btu/hrft-°R 
= length of pipe branch in computer model, ft 
= fluid mass, lb m 

= operator to determine maximum between x and y 
= Nusselt number 
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= fluid pressure, lbf/ft' 

= Prandtl number 

= Prandtl number based on vapor phase properties 
= heat transfer rate at node i, Btu/li 
= radial heat flux, Btu/h-fr 
= tube wall heat flux, Btu/h-ft 2 
= axial heat flux, Btu/h-ft 2 
= tube inner radius, ft 
= Reynolds number 

= liquid- vapor mixture Reynolds number 
= radial direction 
= fluid temperature, °R 
= initial fluid temperature, °R 
= fluid saturation temperature, °R 
= tube wall temperature, °R 
= initial tube wall temperature, °R 
= tune, sec 
= tune increment, sec 
= axial fluid velocity, ft/sec 
= upstream axial fluid velocity, ft/sec 
= mass quality 
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= liquid- vapor mixture correction factor 
= axial direction 

= tube wall characteristic length, ft 
= surface roughness of the tube, ft 

= non-dimensional fluid temperature, (Tf - Tf,o)/(T w ,o- Tf.o) 

= non-dimensional tube wall temperature, (T w - Tf o)/(T w ,o - Tf,o) 
= vapor phase dynamic viscosity, lb„/ft hr 
= fluid density, lb m /ft 3 
= liquid phase density, lbu/ft' 

— density upstream ot a fluid branch in computer model, lb in /ft 
= vapor phase density, lb m /fr 
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= tube wall density, lb n /ft 

= non-dimensional time, (h/p w c w 8)-(t - A cc pfz/m) 

= non-dimensional axial location, nD\\z! riiCf 


Subscripts 

i = computer model node reference 

ij = computer model branch reference 

u = upstream node of computer model branch 


Superscripts 

( I ) = over-dot indicates time rate of change 


Introduction 

The operation of a cryogenic propulsion system requires transfer line chill down 
before establishing a steady flow of cryogenic fluid between various system components. 
Cryogenic transfer line chill down is a transient heat transfer problem that involves rapid 
heat exchange from a solid structure to a fluid with phase change. It is necessary to know 
how long it takes to chill down a given transfer line for satisfactory operation. 

When liquid cryogen, for example hydrogen, at saturation temperature (36.5 °R at 1 
atm) begins flowing through a tube initially at ambient temperature (540 °R) the liquid 
instantly vaporizes near the tube wall. Thus a cross-section ot the flow will have an outer 
vapor ring with a saturated liquid core. As the flow moves downstream, the liquid coie 
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evaporates and the vapor becomes superheated. As the tube wall cools, the liquid coie 
penetrates further and further downstream. Eventually, the tube becomes filled with 
liquid. Reduction in fluid density by vaporization causes the average flow velocity to 
increase signif icantly. Prediction of chill down time requires modeling of these transient 
phenomena and understanding of how they affect heat transfer from the tube wall to the 
flowing cryogen. 

Burke et al. 1 studied chill down of 60 ft., 100 ft. and 175 ft. long, 2 in. 
outer diameter (OD) stainless steel lines by flowing liquid nitrogen. A model to predict 
chill down time was developed by treating the entire line as a single control volume. 

This lumped system provides a simple estimate of chill down time but lacks accuiacy due 
to its broad assumptions and averaging of fluid properties and flow rates over the chill 
down tune. Furthermore, this method cannot be used to calculate instantaneous fluid and 
transfer line wall temperatures. 

Chi 2 looked at chill down of a 26 in. long, 3/16 in. inner diameter (ID) aluminum 
tube using saturated liquid hydrogen as the working fluid. An analytical model of the 
dull down (presented below) was developed under the assumptions of constant flow rate, 
constant heat transfer coefficient, constant fluid properties, homogeneous fluid flow and 
film-boiling dominated heat transfer. The assumption of constant flow rate is not very 
realistic since usually the transfer line inlet and exit pressures are set while the flow rate 
may vary greatly according to the flow condition. The assumptions of constant heat 
transfer coefficient and constant fluid properties are highly restrictive but can provide 
useful estimates of temperature and chill down time. The last two assumptions are 
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idealizations but are widely applicable and produce reasonable results when properly 
applied. 

Steward, et al? modeled drill down numerically using a finite-difference 
formulation of the one-dimensional, unsteady mass, momentum and eneigy equations. 

The model results agree well with experimental results of drill down of a 200 ft. long, 
0.625 in. ID, 0.750 in. OD copper tube using liquid nitrogen and liquid hydrogen as the 
working fluids. Heat transfer coefficients were determined using superposition of single- 
phase forced convection correlations and pool boiling correlations tor both nucleate and 
film boiling. While this type of correlation for film boiling may be acceptable under low 
flow rate/low quality conditions, more recent publications recommend the modified 
Dittus-Boelter type correlation for a wide range of flow conditions. 4 ' 5 ' 6 It is believed that 
improved accuracy and greater generality of application will be achieved tlnough 
implementation of a modified Dittus-Boelter type correlation for flow film boiling. 

Correct modeling of flow film boiling is especially crucial with cryogenic hydrogen since 
film boiling occurs to as low as a tube wall to fluid saturation temperature difference of 
approximately 20 °R resulting in a large portion of the chill down occurring under film 
boiling. 

Complete analytical modeling of the complex heat transfer phenomena occurring 
during chill down is not possible. This paper presents analytical and computer models of 
the cryogenic transfer line chill down process described in five sections. The first section 
is about modeling the chill down process analytically using energy conservation. The 
conservation equation is simplified leading to a closed form solution. This simplified 
case gives insight into parametric effects on the cliill down process as well as serving as a 
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benchmark for the computer solution. The second section discusses theoretical 
development of a heat transfer coefficient tor two-phase heat transfer between the tube 
wall and fluid. The computer model is described in the third section and results from the 
computer model are given in the fourth section. Finally, the fifth section piesents some 
conclusions drawn from the study. 


Analytical Model 

The physical problem modeled is a fluid entering a circular tube at temperature 
T f0 . The tube wall is initially at temperature T w , 0 . Transient heat transfer occurs between 
the tube and the fluid when T w *T f . The flow pattern is approximated as one-dimensional 
with fluid velocity in the axial direction only. Figure 1 shows an energy balance on a 
fluid control volume where the over dots on the energy terms indicate time rate of 
change. 

Conservation of energy tor the control volume states that: 
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Equation (1) can be stated mathematically in cylindrical coordinates neglecting viscous 
dissipation as (see Ref. 7, for example, for derivation) 
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Following the analysis of Chi 2 , Eq. (2) is simplified for application to the problem 
of chill down in cryogenic transfer lines using the following assumptions: 


6 


• Axial conduction in the fluid may be neglected. 

• Flow work may be neglected. 

• Fluid mass flow rate is constant. 

• Heat transfer coefficient (h) is constant. 

• Constant solid and fluid properties 

By expressing q r " in terms of Newton’s law of cooling, and integrating over the 
tube radius, Eq. (2) may be expressed in terms of non-dimensional variable C, [ttDIiz/ 
liicfl as 
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Using a lumped-capacitance model, local wall heat flux may be expiessed in terms of 
non-dimensional variable x [(h/p w c w 5)-(t - A cc ptz/m)] as 
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Equations (3) and (4) together with the initial and boundary conditions 
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Figure 2 shows how non-dimensional tube wall temperature varies with 11011 - 
dimensional time according to a slightly rearranged version of Eq. (5) (see nomenclature 
for definition of 0 W ). For a fixed value of z, T w decreases with time and approaches To 
The effect of varying the parameter C, on the tube wall cooling rate is also shown with the 
four different curves: the tube wall cools more slowly as C, increases. Figure 3 shows how 
non-dimensional fluid temperature varies with non-dimensional tube axial location 
according to Eq. (6). For a fixed value of time, Tf increases with axial location. F 01 a 
fixed value of T f approaches T f , 0 at large time values. 


Heat Transfer Coefficients 

For large differences in tube wall temperature (T w ) and initial fluid temperature 
(Tf, 0 ), film boiling occurs immediately at the entrance of the tube. 9 The mass quality (x) 
of the fluid increases as it travels through the tube eventually becoming superheated 
vapor. In the entrance region of the tube, liquid is surrounded by a vapor film in the so- 
called inverted annular flow pattern. As the liquid-vapor mixture proceeds downstream, 
the liquid phase vaporizes and the flow eventually becomes pure vapor. 

For turbulent single-phase flow, the Dittus-Boelter equation provides a correlation 
of Nusselt number as a function of Reynolds number and Prandtl number 

Nu = 0.023Re°' s Pr° 4 ( 7 ) 
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The above correlation may be applied to a homogeneous model of the two-phase flow 

occurring in the tube using mean fluid properties (see Ref. 10 tor derivation). One such 

••11 1 ^ 

correlation resulting from this type of analysis is that given by Miiopolskn 
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Several researchers have applied a modified Dittus-Boelter formulation to correlate heat 
transfer data. A comprehensive review of correlations of this type may be found in Ref. 
13. 

Since Pr v is close to unity under saturated conditions, this term does not have a 
significant effect on the value of Nu. Therefore, in order to have asymptotic agreement 
of Eq. (8) with Eq. (7) when x = 1 (saturated vapor flow), the exponent on Pr v in Eq. (8) 
is replaced with 0.4. Thus, the proposed liquid- vapor region Nusselt number is 

Nu = 0.023(Re mi J a8 (Pr v )°' 4 (Y) 0 
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Computer Models 

Numerical modeling of boiling heat transfer in a confined tube requires the 
solution of unsteady mass, momentum and energy conservation equations in conjunction 
with a thermodynamic equation of state and correlations for boiling heat tiansfer. The 
Generalized Fluid System Simulation Program (GFSSP) 14 was used to develop a 
computational model of this process. In GFSSP, the conservation equations are first 
expressed in finite volume form for a flow network. A flow network consisting of fluid 
nodes and branches is shown in Fig. 4. At boundary nodes, pressures and temperatures 
are known. At internal nodes and branches, the variables are calculated by solving the 


conservation equations. 

Mass and energy conservation equations are written tor internal nodes to 
determine the pressures and temperatures. The momentum conseivation equations aie 
written for branches to determine the flow rates. The mass conservation equation for a 


typical internal node i can be expressed as: 
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A typical internal node i is connected to the neighboring nodes j through branch lj 
as shown in Fig. 5. 

The energy conservation equation for node i shown in Fig. 5 can be expiessed 
mathematically as shown in Eq. (15). 

1 " ^ 1 r __P_i 

•=^{MAX[-ri'i 1J ,0]t J -MAX[ri') 1J ,0ji 1 }+c] 1 (15) 


m 


i-p- 


pj j 


-in 


V V” A+At 


p j 


At 


j=i 


where 


10 


q,=hA(T„-T„) 


(16) 


and h is determined from Eq. (13). 

The momentum conservation equation tor branch ij can be expressed 


mathematically as shown in Eq. (17). 
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and the Darcy friction factor f is determined from the Colebrook Equation expressed as 
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Equations 14 through 19 are solved numerically. The details of the numerical 
method appear in Ret. 14. All ot the models presented here have oO nodes and use 
hydrogen as the working fluid. Other model parameters are given in Table 1 . 

Equations 14 through 19 are solved numerically in GFSSP. 14 GFSSP uses a 
combination of the successive substitution method and the Newton-Raphson method to 
solve these equations. In this scheme, the mass and momentum conservation equations 
are solved by the Newton-Raphson method while the energy and specie conservation 
equations are solved by the successive substitution method. The underlying principle for 
making such a division is that the equations that are more strongly coupled are solved by 
the Newton-Raphson method. The equations that are not strongly coupled with the other 
set of equations are solved by the successive substitution method. Thus, the computei 
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memory requirement can be significantly reduced while maintaining superior numerical 
convergence characteristics. 

In the iterative loop, mass and momentum conservation equations and the 
equation of state are solved by the Newton-Raphson method. The energy and specie 
conservation equations are then solved by the successive substitution method. Finally, 
the density and other thermodynamic and thermophysical properties and the flow 
resistance coefficient, Kf are calculated. The iterative cycle is terminated when the 
normalized maximum correction, A IIWX , is less than the convergence criterion, C c . A nusx is 

determined from 


=MAX 


y^i 


where Ne is the total number of equations and O is the dependent 
resident mass in node and flow rate in branch). The convergence 
in all models presented in this paper. 


( 20 ) 

variable (pressure and 
criterion is set to 0.01 


Results 

Analytical Model 

Figure 6 shows a comparison of fluid and wall temperatures calculated using the 
analytical model (Eqs. (5) & (6)) with temperatures determined by a GFSSP numerical 
model (Model 1 in Table. 1). As indicated in Table 1, the simulation was done using a 
constant h. Major differences between the analytical and numerical simulations are that 
fluid properties and mass flow rate are allowed to vary in the latter. To approximately 
compensate for variable mass flow rate, a running average of the mass flow rate 
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calculated in the computer solution was iteratively applied to the analytical solution. 
Despite these differences, there is good overall agreement between the two methods. The 
agreement between the analytical and numerical models indicates that the numerical 
solution to the chill down problem is reliable. 

Numerical Model - Superheated Vapor Flow (Model 2) 

To add another level of sophistication towards modeling a fully two-phase flow, 
superheated vapor flow was modeled. In this case, the h for heat tiansfer between the 
tube wall and passing fluid is allowed to vary according to Eq. (13). Since x > 1, Eq. (l->) 
and Eq. (7) are equal for this case. The temperatures vs. time results for this simulation 
are shown in Fig. 7. It is seen that fluid temperature remains relatively low near the inlet 
but increases close to the tube wall temperature near the exit. Chill down is achieved 
after approximately 100 seconds. Figure 8 indicates that h increases with axial location 
until chill down is achieved. Analysis of numerical model output data not shown here 
shows that Reynolds number decreases with axial position but the large increase in fluid 
thermal conductivity with temperature causes h to increase axially, h increases uniformly 
with time near the tube entrance but there is some fluctuation at downstream locations. 
This effect is caused by a combination of increasing Reynolds number and decreasing 
fluid thermal conductivity according to Eq. (7). The effect of Prandtl number on h is 
small. Figure 9 shows that the hydrogen mass flow rate increases until chill down is 
nearly achieved. 

Numerical Model - Subcooled Liquid Flow (Model 3) 

This simulation models the more typical case of cryogenic fluid entering as a 
subcooled liquid during transfer. The subcooled liquid changes to vapor downstream of 
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the tube entrance. The liquid front propagates downstream as the tube cools. Figure 10 
shows temperature vs. time for this simulation. Initially there is significant sensible 
heating of the fluid near the exit of the tube but the temperature quickly drops to 
saturation temperature. Mass quality data (not shown here) from the simulation elucidate 
this trend by showing the fluid to be initially superheated neai the exit and then quickly 
becoming saturated. The fluid is saturated over most of the tube s length duiing chill 
down. The tube wall temperatures display the exact opposite trend in this case compared 
to the single-phase vapor case (Model 2), i.e. the exit cools faster than the entrance. This 
trend is explained by the large increase in h with axial location (see Fig. 11). According 
to Eq. (13), h increases with mass quality. Thus, in the two-phase region, because mass 
quality increases as the fluid moves downstream h also increases axially. Fluctuations of 
h with respect to time are explained with the same reasoning given in the superheated 
case (Model 2). Figure 12 shows that the hydrogen mass flow rate increases and then 
reaches steady-state. Steady mass flow rate for this case occurs earlier on in the chill 
down than in the superheated case when considered as a portion of the total cliill down 
time. 


Conclusions 

Figure 10 indicates that for the case of initially subcooled liquid entering the tube, 
the exit cools faster than the entrance. This phenomenon was also observed by Graham 
et a/. 15 in a series of tube flow film boiling experiments using hydrogen with constant 
wall heat flux. This is opposite the trend for single-phase flow heat transfer. 
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The time to achieve non-dimensional tube wall temperature 0 W = 0 near the tube 
exit (slowest cooling end) is approximately 100 seconds for the superheated hydrogen 
vapor flow case (Model 2). The mass of fluid passed through the tube ovei tliis time 
period is 0.263 lb,,,. In comparison, 2.474 lb„, of fluid is passed through the tube in the 
case of initially subcooled liquid hydrogen over the same time period resulting in © w = 
0.21 near the tube entrance (slowest cooling end). This indicates that the subcooled chill 
down case expends 9.41 times the fluid used in the superheated vapor case to achieve less 
chill down. Tliis result has important implications when reduction in chill down time and 
cryogen expenditure during drill down are primary design considerations. 
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Fig. 1 Energy balance on 


control volume for internal tube flow. 
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Fig. 4 Flow network consisting of nodes and branches. 
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Schematic of fluid nodes, branches and indexing practice 
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Fig. 6 Comparison of GFSSP-calculated wall and fluid temperatures with temperatures 
calculated by analytical solution (Model 1: superheated hydrogen, eonstant h). 
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Fig. 7 Temperature vs. time at different axial locations (Model 2: superheated hydrogen 
variable h). 
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Fig . 8 Heat transfer coefficient vs. iime a. differen, axial locations (Mode. 2: superheated 
hydrogen). 
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Fig. 10 Temperature vs. time at different axial locations (Model 3. 
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hydrogen, variable h). 
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Fig. 12 Mass flow rate vs. time (Model 3: initially subcooled hydrogen). 




